package edu.gmu.atelier;

/**
 * Performs QR factorization on the passed matrix using Householder
 * reflections.
 * 
 * @author James H. Pope
 */
public class QRHouseholderFactor implements QRFactor
{
    private Matrix a  = null;
    private Matrix q  = null;
    private Matrix r  = null;
    
    public QRHouseholderFactor( Matrix a )
    {
        this.a = a;
        
        /*
         * Go look at http://www.math.iit.edu/~fass/477_577_handouts.html
         * Chapter 4 seems to provide a good explanation of QR Factor.
         * 
         * Also MIT
         * http://persson.berkeley.edu/18.335/
         * Specifically
         * http://persson.berkeley.edu/18.335/lec6handout6pp.pdf
         * and
         * http://people.maths.ox.ac.uk/trefethen/text.html
         */
        this.decompose(a);
    }
    
    public void decompose( Matrix a )
    {
        boolean verbose = false;
        
        this.r = a.copy();
        
        /*
         * Using House Holder reflectors
         */
        int n = r.getCols();
        int m = r.getRows();
        
        // Remember the u's, they become the matrix Q
        this.q = r.getFactory().makeIdentity(m); // presumably the ones are overwritten?
        
        // Determine the maximum number of iterations required
        // Note that this allows m=n, m<n and m>n matrices
        int t = Math.min( m-1, n );
        if(verbose)System.out.println("\n\nDecomposition with "+t+" iterations.");
        //int t = 1;
        //int t = Math.min( m-1, n );
        Matrix i = r.getFactory().makeIdentity(m);
        //System.out.println(""+m+" x " + m + " = "+i);
        for (int k = 0; k < t; k++)
        {
            //Matrix minor = r.copy(k, t-k);
            Matrix minor = r.copyMinor(k, k);
            if(verbose)System.out.println("Minor="+minor);
            
            Vector x = minor.getVector(0);
            
            Element kth = x.getElement(0);
            Element norm = x.norm2().sqrt(); // this is definitely positive
//            if( kth.isNegative() )
//            {
//                norm = norm.negate();
//            }
//            //norm = norm.negate();
            if(verbose)System.out.println("||x||="+norm);
            
            // Ugly, find better way, make e1 same length as x
            Matrix fi = r.getFactory().makeIdentity(minor.getRows(), minor.getCols());
            //System.out.println("Fi="+fi);
            Vector e1 = fi.getVector(0);
            
            /*
             * NOTE: This is not complete for Complex Matricies, need to test
             * If real w always equals 1, therefore will be 2 times the vvt
             * If complex, w = x^T.u / u^T.x
             * I believe this supports both?
             */
            Vector xe1 = e1.mult(norm);
            if(verbose)System.out.println("xe1="+xe1);
            Vector v = x.add( xe1 ); // can be add or subtract?
            //this.q.getVector(k).copyElements(k, u); // nor sure correct?
            if(verbose)System.out.println("v="+v);
            Vector u = v.normalize();
            if(verbose)System.out.println("u="+u);
            Element one = kth.getOne();
            Element vtx = u.dot(x);
            Element xtv = x.dot(u);
            Element w = xtv.mult( vtx.inverse() );
            Element multiplier = one.add(w);
            if(verbose)System.out.println("Multiplier="+multiplier);
            Matrix twovvt = fi.mult(u, u).mult(multiplier);
            
            Matrix f = fi.sub(twovvt);
            if(verbose)System.out.println("F="+f);
            //System.out.println(" I="+i);
            Matrix qk = i.copyElements(k, k, f);
            //System.out.println("Qk="+qk);
            q = q.mult(qk);
            
            // First part of Q is the I identity matrix
            if(verbose)System.out.println("Qk="+qk);
            this.r = qk.mult(r);
            if(verbose)System.out.println("\nIteration: "+k+" R="+r);
        }
        //System.out.println("Q="+q);
        //System.out.println("R="+r);
    }
    
    public Matrix getA()
    {
        return this.a;
    }
    
    public Matrix getQ()
    {
        return q;
    }
    
    public Matrix getR()
    {
        return r;
    }
}
